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According to the Jarzynski theorem, equilibrium free energy differences can be calculated from 
the statistics of work carried out during non-equilibrium transformations. Although exact, this 
approach can be plagued by large statistical errors, particularly for systems driven strongly away 
from equilibrium. Recently, several approaches have been suggested to reduce these errors. In this 
paper we study the efficiency of these methods using two models for which analytical solutions exist. 



I. INTRODUCTION 

The calculation of free energies is central to many ap- 
plications of computer simulations ranging from the cal- 
culation of ligand affinities to the study of phase equi- 
libria in condensed materials. Since the free energy is a 
quantity related to the phase space volume available to a 
system, its calculation is non-trivial in most interesting 
cases and in the past decades considerable effort has been 
devoted to developing efficient free energy computation 
techniques [2,0]. 

Recently, Jarzynski has shown how free energies can be 
determined from non-equilibrium simulations [3( • In this 
so called fast switching or fast growth method a control 
parameter coupled to the Hamiltonian of the system is 
switched from an initial to a final value at a finite rate. 
By changing the external parameter, the work W is per- 
formed on the system. As a consequence of the second 
law of thermodynamics the average of this work is larger 
than the free energy difference AF between the equilib- 
rium states corresponding to the final and initial value of 
the control parameter: 



(W) > AF. 



(1) 



Here, the angular brackets denote an average over path- 
ways starting from the equilibrium state corresponding 
to the initial value of the control parameter. The equal 
sign in this inequality, known as Clausius inequality or 
maximum work theorem, holds only if the external pa- 
rameter is switched reversibly. 

In 1997 Jarzynski demonstrated that the Clausius in- 
equality can be turned into an equality by averaging the 
work exponential rather than the work itself [H, Q : 



-pw\ 
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(2) 



where (3 = \/k^T reciprocal temperature. The Jarzynski 
equation ^ has been proven for several types of dynam- 
ics ranging from Newtonian and thermostatted dynamics 
to Langevin and Monte Carlo dynamics [1, 0, H,0, S Hi- 
Essentially, the theorem is valid if the dynamics pre- 
serves the canonical distribution for fixed control param- 
eter 0,3. 

Remarkably, the validity of the Jarzynski equality does 
not depend on how the control parameter is switched 



from its initial to its final value. (In particular, it is not 
even necessary to maintain a fixed protocol throughout 
the simulation or the experiment as long as the initial and 
final value of the control parameter are fixed.) Thus, the 
switching process can be carried out at arbitrary speed. 
For infinitely slow switching, Jarzynski's method reduces 
to Kirkwood's coupling parameter method, or thermo- 
dynamic integration [9j], while for infinitely fast switch- 
ing one obtains Zwanzig's thermodynamic perturbation 
method (To| . Thus, these two well known free energy cal- 
culation methods can be viewed as limiting cases of the 
fast switching procedure. 

Whether the freedom derived from the arbitrariness 
of the switching protocol can be used to develop more 
efficient free energy calculation methods is an interest- 
ing and practical question. A recent analysis indicates 
that the statistical errors encountered in evaluating ex- 
ponential averages such as that of Eq. @ prevent fast 
switching methods from exceeding the efficiency of con- 
ventional methods such as umbrella sampling or ther- 
modynamic integration [ill ]. This conclusion seems to 
remain true even if path sampling techniques are used to 
focus on those non-equilibrium pathways that yield the 
larg est contribution to Jarzynski's exponential average 
[l2t [l3|. Note that these path sampling approaches can 
also be combined with the large time step method for 
fast switching simulations fbij . 

In this paper we analyze the efficiency of several 
fast switching methods including the path sampling ap- 
proaches mentioned above for two analytically solvable 
models: a particle dragged through a viscous fluid in a 
harmonic trap [Hj] and an ideal gas in an expanding pis- 
ton [lj| . For these two models the work distributions can 
be calculated analytically for arbitrary switching speeds. 
Using these work distributions, we determine the statis- 
tical errors and, from them, the computational efficiency 
that one would obtain in fast switching simulations of 
these systems. The results presented in this paper con- 
firm our previous conclusions based on actual simulations 
[ril ]. Namely, fast switching simulations do not outper- 
form umbrella sampling and thermodynamic integration 
in most interesting cases. 

The remainder of the paper is organized as follows. 
Simulation methods are described in Sec. [TTJ Expressions 
for error and efficiency estimates are presented in Sees. 
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IIIII and HVl In Sec. [V] these expressions are then applied 
to the two model systems. Conclusionss are provided in 

sec. rvn 



II. FAST SWITCHING METHODS 

In this section we review the fast switching methods 
whose efficiency we will analyze in later sections. To 
set the notation, consider a system with a Hamiltonian 
Tt(x, A) that is a function of the state of the system x 
and the external control parameter A. Depending on the 
system, x may include positions and momenta of all par- 
ticles or the positions only. The free energy of the system 
at temperature T and fixed control parameter A is given 
by 



F x (P) = -k B TlnQx((3), 
where the normalizing factor 



Qx(J3) 



dx 



-0H{x,X) 



(3) 



(4) 



is the partition function of the system (up to a constant 
factor). The free energy difference between two states A 
and B corresponding to two different values of the control 
parameter, A^ and As, respectively, is given by 



AF = F\ B - F\ 



-fceT In 



Q\ B 



(5) 



State A is described by the Hamiltonian Ti.(x, Xa) and 
state B by H(x,Xb)- We now imagine that the con- 
trol parameter A is continuously switched from the initial 
value A^ to the final value A^ according to some proto- 
col A(i) within a total time r. As the control parameter 
changes, the system evolves in time tracing out a tra- 
jectory x(t). Here, x(t) denotes a complete trajectory 
describing the system from time t = to time t = r. 
The change in energy that is due to the changing control 
parameter is the work W carried out on the system along 
the particular trajectory x(t): 



W[x(t),X(t)]=£dtX(t)^- 



(6) 



=x(t) 



where X(t) — dX(t)/dt. This notation emphasizes that 
the work W depends both on the specific trajectory fol- 
lowed by the system as well as on the particular protocol 
used to switch the control parameter. In path integral 
notation, the Jarzynski equality can be written as: 



-/3AF 



Vx(t) P[x(t),A(i)]e-^ w[x(i) ' A(i)1 . (7) 



Here, J T>x(t)... denotes a summation over all trajecto- 
ries and P[x(t), X(t)} is the probability density of tra- 
jectory x(t), along which the work W[x(t), X(t)] is per- 
formed during the switching process. The probability 



density V[x(t), X(t)] includes the canonical distribution 
p(xo) = exp[— (3H(xq, Xa)]/Qx a of the initial conditions 
xq. For deterministic dynamics, the trajectory x(t) and 
the work W[x(t) 7 X(t)} are fully described by the initial 
point in phase space xq. In this case, the integral over 
all trajectories can be written as an integral over phase 
space rather than an integral over trajectory space. In 
the following, we will often omit the explicit dependence 
of W^a^i)] and V[x(t)] on X(t) to simplify the notation. 



A. Straightforward fast switching 

The Jarzynski equation (|7|) justifies the following algo- 
rithm to estimate the free energy difference between the 
two equilibrium states A and B corresponding to Xa and 
Xb, respectively. First, one generates N initial conditions 
canonically distributed with respect to H.(xo, Xa)- This 
can be done with Monte Carlo sampling or constant tem- 
perature molecular dynamics. From each of the N initial 
condition Xq one then generates a trajectory of length r 
by integrating the appropriate equations of motion and 
calculates the work performed on the system along 
that trajectory. From this sample of N trajectories an 
estimate AFn of the free energy is then determined: 



AF 



iV 



-k^T In 



1 N 

-Y 



(8) 



i=i 



As the trajectory sample is finite, this free energy esti- 
mate contains errors [171 ] that are discussed in detail in 
Sec. IIIII Due to the highly non-linear behavior of the 
exponential function these errors may be severe, often 
precluding an accurate free energy calculation particu- 
larly in the case of large switching rates ll| . The reason 
for these large deviations is best perceived by considering 
the work distribution P(W): 

P(W) = Jvx(t)V[x(t)]S(W -W[x(t)}). (9) 

In terms of P(W), the Jarzynski equation can be rewrit- 
ten as an integral over the work: 



= -/3AF 



dWP(W)e 



-pw 



(10) 



In a fast switching simulation as described above, most 
trajectories have work values near the values for which 
P(W) is a maximum. For large switching rates, P(W) 
may be centered at large work values, such that for 
typical work values exp(— j3W) is a very small number 
yielding a vanishing contribution to the integral in Eq. 
(flTj|) . The work values leading to the important contri- 
butions to the integral, on the other hand, may occur 
very rarely [l8j . This situation, familiar from the calcu- 
lation of chemic al p otentials with Widom's particle in- 
sertion method [19( or from Zwanzig's perturbation ap- 
proach [lCj| . may lead to extremely large statistical in- 
accuracies. These statistical difficulties can easily offset 
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the advantage of using short, computationally inexpen- 
sive non-equilibrium trajectories. 



B. Work biased thermodynamic integration 

To avoid large statistical errors occurring in the 
straightforward application of the Jarzynski relation, Sun 
proposed a method based on a thermodynamic integra- 
tion in trajectory space [l2|. The basic idea of this 
method is to introduce a bias that favors the sampling 
of the rare but important trajectories that mostly con- 
tribute to the exponential work average. This is done 
by introducing a free energy AF(a) that depends on the 
parameter a: 

e -/JA#(a) = fvx(t)V[x(t)]e- paW ^. (11) 



For a = 0, the parameter dependent free energy vanishes 
because P[x(i)] is normalized, AF(0) = 0. For a = 1, on 
the other hand, the parameter dependent free energy is 
identical to the original free energy, AF(1) = AF. One 
can therefore calculate AF by taking the derivative of 
AF(a) and integrating it from to 1: 



f 1 , dAF(a) 
AF = da- 
Jo 



da 



(12) 



To carry out this procedure, which is the path space ver- 
sion of the thermodynamic integration method one 
needs the derivative of AF(a) with respect to a: 



dAF(a) _ JVx(t)V[x(t)]e-^ aW ^W[x(t)} 



da 



JVx{t)V[x(t)]e-^ w ^ 



(13) 



The right hand side of this equation can be viewed as 
the average (W) a of VF[a;(i)] over the work biased path 
ensemble 



V a [x(t)]=V[x(t)}e-^ w ^/Z a , 



where 



Z a = J T>x(t)V[x(t)]e 



PaW[x(t)] 



(14) 



(15) 



In terms of the work average (W) a , Eq.fl2|) can be 
rewritten as: 



AF 



da (W) c 



(16) 



To calculate the free energy difference AF using this 
equation it is first necessary to determine (W) a for dif- 
ferent values of the parameter a between and 1. Then, 
the integral on the right hand side is determined numer- 
ically from these values. As shown by Sun [l2T ]. path 
sampling techniques, orig inally developed to study rare 
but important events 20 , [2l[ , can be used to sample the 
work biased ensemble and calculate the averages (W) a . 



To understand why the work biased thermodynamic 
integration yields accurate free energy estimates, it is in- 
structive to consider the work distributions P a (W) in the 
work biased ensembles: 



(17) 



P a (W) = I Vx(t)V a [x(t)}6(W -W[x(t)}) 

P(W) exp(-0aW) 
J dWP{W) aip(-f3aW) 



For a = the parameter dependent work distribution 
equals the original work distribution: 



P (W) = P(W). 



(18) 



In the other limit, at a = 1, the parameter dependent 
work distribution is proportional to the integrand in Eq. 
(EQD: 



P^W) = P{W)e 



-@{W-AF) 



(19) 



For values of the parameter a between and 1 the work 
distribution in the work biased ensemble is intermediate 
between these two limiting cases. Thus, both the typical 
and the dominant work values are generated with com- 
parable likelihood as a is changed from to 1 yielding a 
good statistical accuracy. 



C. Work biased umbrella sampling 

An alternative way to enhance the sampling of rare but 
important work values was recently proposed by Ytre- 
berg and Zuckerman [l3| and Athenes [23]. The basic 
idea of this approach, called single-ensemble nonequilib- 
rium path-sampling by Ytreberg and Zuckerman, is to 
introduce an explicit biasing junction (or umbrella func- 
tion) 7r[a;(t)] depending on the path x(t). With this um- 
brella function the Jarzynski equality ([2]) can be rewrit- 
ten as: 



-f)AF 



JVx(t)V[x(t)]n[x(t)] 




-«»[-W] " 




w[x(t)] 




JVx(t)V[x(tMx(t)} 


1 

_ir[x(t)]_ 





(e-^W'l/iKt)]), 



(20) 



Here, (• 
semble 



(l/7T[x(t)]) v 

)n denotes an average over the biased path en 



V{x(t)Mx(t)} 
JVx(t)V[x(t)]n[x(t)} 



(21) 



For a path observable A[x(i)] depending on the pathway 
x(t) the average in the biased path ensemble is given by: 



{A[x(t)])* 



JVx(t)V[x(t)]n[x(t)]A[x(t)] 
fVx(t)V[x(t)]*[x(t)] 



(22) 



4 



Such averages can be calculated using path sampling 
techniques [H, HO, HH, HH ■ Since the umbrella function 
is introduced to guide the sampling to the relevant work 
values, it suffices to let 7r[a;(f)] depend on the work only, 
7r[a;(f)] = 7r[Vy(a;(i))]. In this case, the biased average 
of a function A[W(x(£))] that also depends on the work 
only can be written as 



_ JdW P(W)ir(W)A(W) 
[ U ~ J dW P(W)n(W) 



(23) 



This expression will be useful in the next section, in which 
we discuss the statistical errors occurring in fast switch- 
ing simulations. 

In a work biased umbrella sampling simulation the 
choice of the umbrella function 7r[x(t)] is crucial. To 
obtain a good accuracy in the free energy estimate, it 
is necessary that the bias function has a good overlap 
with the unbiased work distribution P(W) as well as 
with the integrand of Eq. l[T0"|). P(W) exp(-(3W). Ytrc- 
berg and Zuckerman found that the umbrella function 
ir(W) = exp(— /3W/2) leads to considerable efficiency in- 
creases compared to the results of straightforward fast 
switching simulations. Another possibility would be to 
carry out umbrella sampling simulations with partially 
overlapping windows [l[ . 

Note that other non-Boltzmann sampling techniques 
such as flat histogram sampling 12311 . multicanonical sam- 
pling [U, or parallel tempering [25| may also be used to 
improve the convergence of the exponential average of 
the Jarzynski equation. These techniques, however, are 
not considered in the present article. 



III. ERROR ANALYSIS 

In this section we review the expressions required to 
estimate the statistical errors arising in fast switching 
simulations. While straightforward fast switching and 
work biased umbrella sampling can be treated in the same 
way, separate expressions are required for the work biased 
thermodynamic integration scheme. 



A. Straightforward fast switching and work biased 
umbrella sampling 

To set the notation we reproduce here the expressions 
derived in Ref. ll|. All expressions are valid both for 



work biased umbrella sampling as well as for straightfor- 
ward fast switching. The latter case follows by setting 
the umbrella function equal to unity, 7r[a;(i)] = 1. 
For convenience we define 



X[x(t)} 



-PW[x(t)] 



and Y[x(t)} 



1 



(24) 



Tr[x[t)] ir[x[t)Y 
The free energy estimated according to Eq. (|2"0|) from N 



pathways is 



AF 



-fceTln 



X 
V 



N 



(25) 



JV 



where Xn and Y n are averages over the N trajectories 
of the simulation. 



X^lf>(*) and F^lf>«. 



(26) 



; = 1 



Here, X« and Y« arc the values of X and Y associated 
with the i-th path sampled from the biased ensemble, 
P*[x(t)]. 

Due to the non-linearity of the logarithm, this estima- 
tor for the free energy is biased, i.e., the average over 

many realizations AF N of the free energy estimator, 



(AF 



M 

A/-00 M N 
3=1 



(27) 



differs from the true free energy difference. This devia- 
tion, 



b N = (AF N ) - AF, 



(28) 



is called the bias. In addition to this systematic devi- 
ation, the estimator AFn is also affected by statistical 
errors which are quantified by the variance, 



'N 



\AF 



N 



(AF 



((AF 



(AF N ) 2 . (29) 



The total mean squared deviation of the estimator AFn 
from the true free energy difference AF has contributions 
from the bias and the statistical error: 



~N 



([AF N - AF} 2 ) = b 



2 
N 



(30) 



We now consider the limit of a large sample size N and 
first calculate the bias b]y. For large N, the deviations 
SXn and SYm from their averages Xn and Yjv, 



SX 



N 



X 



N 



(X)tt and 5Y N =Y N -(Y)„, (31) 



are small compared to (X) n and (Y) n , respectively. Ex- 
pansion of the logarithm in Eq. (|25| into a Taylor series 
around (X) / (Y) and truncation after the quadratic term 
(the linear terms vanish after averaging) yields the bias 



>N 



2N 



((sxy< 



((SY) 2 ), 



(32) 



where the fluctuations SX and SY are given by 

SX = X - (X) n and 5Y = Y - (Y) n . (33) 



Interestingly, Eq. (|32|) implies that the bias can be made 
to vanish by choosing an appropriate umbrella function. 

For the variance ct/v and the mean squared error 
the Taylor series of the logarithm in Eq. ([23)) can be 
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truncated after the linear term. In this approximation 
the bias vanishes and the variance and the mean squared 
error are identical: 



>N 



N 

TO 



((5YY 



(Y)l 



- 2 



(6X6Y)„ 



.(34) 



In the absence of an umbrella function tt (or. in other 
words, for tt = 1) these expressions describe the errors 
encountered in straightforward fast switching simulations 
and are identical to the ones derived in Refs. [l7| and 

M- 

All terms occurring in the above expressions for the 
bias the variance ctjv, and the mean squared error 
can be calculated as integrals involving the work distri- 
bution P(W): 



I = J dW P(W)n(W), 



1 

7 
i 

T 
i 

7 

i 

7 

i 

7 



(X) n = 9 / dW P{W) cxp(-pW) 
(XY) n 



dW P(W) 
dW P{W) 



exp(-2/?W) 
w(W) 



dW P(W) 



ir(W) 
exp(-/?W) 



ir(W) 



(35) 
(36) 
(37) 
(38) 

(39) 
(40) 



In section IVl we will use these equations to calculate the 
expected errors for two systems with analytically known 
work distributions. 

In the following we will also consider the errors in the 
free energy if the process is carried out in reverse di- 
rection, i.e., if the system starts from equilibrium initial 
conditions corresponding to the final value of the control 
parameter and the system evolves under a time-inverted 
switching protocol. In general, these errors differ from 
those of the forward process. However, in the case of 
umbrella sampling with bias tt(W) = exp(— (3W/2) the 
errors of the forward and reverse process are identical. 
To show this we consider the work distribution P R {W) 
of the reverse process, which is related to the work dis- 
tribution of the forward process by the Crooks identity 

sua 

-0W+0AF 



Pr{-W) = P{W)e 
Insertion of Eq. (gTJ) into Eqs. fl35]) to (gDJ) yields 



(X 2 



(Y 2 



(X)l. R 



(Y)l F 



(X 2 



'ir,F 



(Y)ln 
(XY) n , R 



(X)l, F ' 
(XY)^f 



(X) ViR (Y)„ iR (X) 7r . F (Y) nt 



(41) 

(42) 
(43) 
(44) 



Here, {■■} 7 :,f and (-)^,r denote averages with umbrella 
function n(W) for the forward and the reverse process, 
respectively. Thus, the fluctuations of X in the forward 
process are identical to the fluctuations of Y in the re- 
verse process and vice versa. The correlation term is the 
same for both directions. Inserting these results into Eq. 
(|34"1) one finally obtains 



-N,F 



(45) 



i.e., the mean squared errors for the forward and reverse 
process are identical. In the same way one can show that 
also for the 1/P-bias the errors are the same in forward 
and reverse direction. 

For straightforward fast switching Eq. (|34|) is very sim- 
ilar to an expression obtained recently by Jarzynski that 
relates the number of trajectories required to obtain ac- 
ceptable accuracy to the dissipative work of the reverse 
process [l8[ . To establish this similarity we use Eq. |34|) 
to determine the number N^t of trajectories required to 
obtain a statistical error of thermal magnitude k^T: 



J dWP{W) exp(-2/3W0 
exp(-2/3AF) 



1. 



(46) 



Using the Crooks identity, Eq. (|46|) can be written as 
/ dWP R (-W) exp(-/W) 



kT 



exp(-AF) 



(47) 



where we have neglected the term -1 on the right hand 
side of Eq. (j46|) (in any interesting case N is much larger 
than 1). Changing the integration variable from W to 
— W one obtains: 



kT 



dWP R (W) exp[/3{W + AF)) 



(expG0W£)) fl , 



(48) 



where W R = W + AF is the dissipative work for the re- 
verse transformation and (• • ■) R indicates an average over 
realizations of the reverse process. This equation implies 
that the number of fast switching trajectories required for 
an accurate free energy estimation depends on the work 
dissipated during the reverse process, as noted earlier by 
Jarzynski [l8[ . Remarkably, for highest accuracy of the 
free energy estimate the process should be carried out in 
the direction that yields the larger exponential average of 
the dissipative work [l8| . Of course, the same conclusion 
can be drawn from 



NkT = (cxp(-2/3W d )), 



(49) 



which also simply follows from Eq. (|4"6"|) and where W d = 
W — AF is the dissipative work for the forward process. 

To simplify Eq. (|4"8|) , one is tempted to replace the 
average of the exponential with the exponential of the 
average (l8j . 



(exp(f3W d )) R ^exp((3W R ), 



(50) 
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where W R = J dWP R (W)(W + AF) is the average dissi- 
pative work for the reverse process. This approximation, 
however, is not valid in general. As shown for an illus- 
trative example in Sec. IV A| it is likely to underestimate 
the required number of trajectories that may be off by 
orders of magnitude from the true value. 



B. Work biased thermodynamic integration 



According to Eq. p^|) . Sun's work biased thermody- 
namic integration method requires the calculation of an 
integral over the parameter a. To determine this inte- 
gral, the average (W) a is first calculated for L discrete 
values otj of the parameter a between and 1 in L sep- 
arate simulations. Then, the integration is carried out 
numerically using for instance the trapezoidal or Simp- 
son's rule. Imagine now that the average work (W) a is 
approximated by averaging over M trajectories at each 
value of the parameter a: 



(W) a « W M {a) 



1 M 







(51) 



i=0 



Here, the argument a in W m{o) indicates that the tra- 
jectories over which the average is carried out are sam- 
pled from the trajectory ensemble corresponding to the 
parameter a. Using the trapezoidal integration rule the 
estimate for the free energy is then given by 



AF 



N 



1 

3=1 



W M {0Lj), 



(52) 



where N = ML is the total number of trajectories gener- 
ated in the simulation. Neglecting errors in the numerical 
integration over a this free energy estimate is bias free: 



(AF N ) =jYl(Wu(a j )) = jf^(W) aj =AF, 



(53) 



3 = 1 



3=1 



where we have used that (Wm^j)) = (W) aj . In the 
above equation the angular brackets (■ • •) denote an av- 
erage over many independent realizations of the thermo- 
dynamic integration procedure. 

To determine the variance afj of the free energy esti- 
mate, which in this case is equal to the mean squared 
error t 2 N , we note that the variance of a sum of random 
variables equals the sum of the variances of the individ- 
ual variables. Accordingly, the variance of AFn is given 
by 



1 

I 2 



L 

£* 

3=1 



IV 



where 



w 



(a,) = ([W M (a 3 ) ~ (W) aj } 2 ) 



(54) 



(55) 



is the variance of the average W m(c*j) at parameter value 
Oij. For statistically independent work values this vari- 



ance is given by 



w 



i a i) 



M 



[(W 2 



(56) 



Rewriting the sum over j in Eq. 
a we finally obtain: 



' JV 



1 

N 



54]) as an integral over 



(57) 



[(W 2 ) a - (W)%] da. 



Note that here we have assumed that at each value of 
the parameter a, the same number M of trajectories are 
generated. To minimize the variance a 2 ^, however, it is 
advantageous to sample more trajectories at values of a 
where the work variance is large and less where it is small. 
For simplicity we do not consider this case here. 

For a given work distribution P(W) the work variance 
in the work biased ensemble can be calculated from 



and 



(W) a 



(W 2 ) a = 



J dWP{W)e~ al3W W 
J dWP{W)e~ 



a/3W 



J dWP{W)e- al5W W 2 
j dWP{W)e- a f )w 



(58) 



(59) 



By numerical integration of Eq. ([57)) one can then calcu- 
late the variance cr^ of the free energy estimate AFn- 

Interestingly, also the error of the thermodynamic in- 
tegration method for the reverse process is identical to 
that for the forward process. The average of an arbitrary 
function A(VF) for a particular a of the reverse process 
is 

(A(m) _ JdWP R (W)e-^ w A(W) 



/ dWP R (W)e- a 



Inserting the Crooks identity, Eq. (|4"Tj) , into Eq. (f6T)]) we 
obtain 

. . f dWP{W)e-^ 1 -^ w 'A(-W) , N 

= (A(-W))i- a ,p. (62) 

Thus, the work variance for the forward process at a 
is equal to the work variance for the reverse process at 
(1-a): 

(W 2 ) a , R - (W)l >R = (VF 2 )i- q ,f - (W) 2 _ a;F (63) 

Integration over a according to Eq. (|5T[) with a change 
of variables from a to 1 — a in the reverse case then 
demonstrates that the statistical errors for the forward 
and reverse process are indeed equal: 



N,F 



N,R- 



(64) 
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IV. EFFICIENCY 

While the error obtained from a simulation with N 
trajectories is an interesting quantity, the criterion most 
relevant for the practitioner is the computational cost 
required to achieve a given accuracy in the free energy. 
For the estimation of this cost one needs to take into 
account that long trajectories are computationally more 
expensive than short ones. In order to do that we note 
that for all algorithms studied here the error in the free 
energy can be written as: 



4 



K 



(65) 



where k is a constant that does not depend on the number 
of trajectories used in the calculation. Hence the number 
NkT of trajectories required to obtained an accuracy of 
k^T in the free energy is given by 



N kT = 



klT>- 



(66) 



Here, wc have set the target error to a value of k^T 
because variations of the free energy much smaller than 
kftT are irrelevant in most cases. 

Since the computational cost of a trajectory is pro- 
portional to its length r a meaningful definition of the 
computational cost is given by 



Ccpu = NkTT 



fc|T* 



(67) 



Ccpu is the computational cost of a free energy calcula- 
tion with accuracy k^T measured in units of the compu- 
tational cost to generate a trajectory of length r = 1. 
While for straightforward fast switching and umbrella 
sampling the computational cost is 



CPU 



{{SXf), , ((<5F) 2 ). „ (8XSY)« 



(Y)l 



(68) 

for the work biased thermodynamic integration scheme 
is given by: 



Ccpu 



P 2 r 



(W)l)da. 



(69) 



V. RESULTS 



A. Particle pulled through a viscous fluid 

1. Model and work distribution 

The first model we consider is a particle in one dimen- 
sion pulled through a viscous liquid at temperature T by 
a harmonic trap with force constant k moving at constant 
speed v for a certain distance L. The translation of the 
trap by the distance L takes a time r = L/v. The motion 
of the particle is modelled by a Langevin equation in the 
overdamped limit with friction constant 7: 



djj k . . 

_ = --(*-««)+„, 



(70) 



where x is the position of the particle and r\ is a Gaus- 
sian random variable uncorrelated in time, (r]{t)r](t')} = 
(2ksT /^i)8(t — t'). For this model, Mazonka and Jarzyn- 
ski have computed the work distribution analytically 
finding a Gaussian distribution (Til ]: 



P(W) = 



1 



2tt(Tw 



exp 



(W - Wf 



with mean 



W = jLv 



and variance 



'W 



} jLv 



1 



v 
~k~L 



-kL/wy 



1 



w. 



(71) 



(72) 



(73) 



Since the free energy of the system does not depend on 
the position of the harmonic trap, AF = for this trans- 
formation. The average work, however, is larger than 
zero due to the energy dissipated into the bath and only 
in the limit of infinitely slow switching does the average 
work vanish. For slow switching, i.e., for v kL/j or 
r ^> 7/fc, the average work is proportional to the pulling 
velocity, W = jLv, as predicted by linear response the- 
ory. In the limit of infinitely fast switching, i.e. for 
v 3> or t <C 7/fc, the average work converges to 

a constant value, lim^oo W = kL /2. Accordingly, the 
variance erj^ = 2k^T^Lv for slow switching. In the fast 
switching regime the variance reaches a constant value, 
°~w = k^TkL 2 (see Fig. [T]). All results shown in this 
section were obtained for k^T = 1, 7=l,fc = l and 
L = 5. 



If the work distribution P(W) is known for a specific 
switching process, the equations of the previous sections 
can be used to determine the expected accuracy of the 
free energy calculated with the methods of Sec. |TTJ Here 
we do that for two models, for which the work distri- 
butions are known analytically. In particular, we study 
the computational efficiency of these methods in terms of 
the computing time required to obtain a given accuracy 
in the calculated free energy. 



2. Efficiency 

For this model we have calculated the expected errors 
for straightforward fast switching (i.e., ni(W) = 1) as 
well as for the two bias functions ^{W) = exp(— (3W/2), 
and tt 3 (W) = 1/P(W). The bias function n 2 (W) = 
cxp(— (5W/2) was suggested by Ytreberg and Zuckerman 
and is designed to improve the sampling of low work val- 
ues [lH . Some attention is needed when considering the 
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FIG. 1: Work variance a^, for a particle pulled through a 
viscous fluid by a harmonic trap moving at constant speed v. 
The average work is given by W = /3a^/2. 



umbrella function 773 (W). Since this function biases the 
simulation in a way to make the work distribution flat, 
tts(W) can have the form 1/P(W) only in a range with 
finite size. We therefore define 

f l/P(W min ) for W < W min , 
tt 3 (W) = { 1/P(W) for W min <W< W max , 
{ l/P(W max ) for W > W max . 

(74) 

The limits for this range, W m i n and W max , need to be 
selected such that all important work values are included 
in the sampling. Of course, in the general case the work 
distribution P(W) required to implement the umbrella 
function ir 3 (W) = 1/P(W) is not known in advance. 
However, methods such as flat histogram sampling [23[ 
can be used to determine this bias itcratively. Here, we 
study this case in order to evaluate the efficiency of the 
work biased umbrella sampling approach if a good um- 
brella function is available. 

For the Gaussian work distribution of the pulled har- 
monic oscillator and the bias functions considered here 
all integrals in Eqs. |35|) to (|40|) can be computed ana- 
lytically. These analytical expressions can then be used 
to calculate the expected errors for straightforward fast 
switching, 



N 



(f<& - l) 



(75) 



work biased umbrella sampling with the exponential um- 
brella function tt(W) = exp(— f3W/2), 



N 



and with the umbrella function n 3 (W) = 1/P(W), 



2 _ k 2 B T 2 U 



N y/Trcr w 



1 _ e -/3 2 ^-/4 



(76) 



(77) 



the interval U in which the work distribution is required 
to be flat we use the distance between the maximum of 
the work distribution P(W) and the maximum of the 
Jarzynski integrand P(W) exp(— j3W). For the Gaussian 
work distributions considered here this distance is /3cr^ 
and we set U — (icily accordingly. With this choice, the 
statistical error of the \ jP bias becomes: 



k\T 2 



flaw 



N 



1 



(78) 



The expected errors are shown in Fig. Has a function 
of the switching rate v. Since in all expressions for the 
error the number TV of trajectories factors out in a simple 
way, we have represented k 2 = t 2 N N as a function of v 
in this figure. The number N^x of trajectories required 
to obtaine an error of k^T is related to k 2 simply by 
N kT = n 2 /k B T. 

We also calculated the expected errors for Sun's work 
biased thermodynamic integration scheme. Direct eval- 
uation of the work distribution P a (W) [see Eq. (|17p] 
in the work biased ensemble shows that for a Gaussian 
distribution P(W) the variance of the work biased dis- 
tribution P a (W) is the same for all values of a. In this 
case Eq. (|57| yields 



1 

N 



((w 2 ) - (wo 2 ) 



1 

N 



(79) 



While the variance in the work biased ensemble does not 
change as a function of a, the average work (W) a does: 



(W) a = (W) 



(80) 



The expected errors for the work biased thermodynamic 
integration method are also shown in Fig. [5] as a func- 
tion of the switching rate v. While for straightforward 
fast switching and for umbrella sampling with exponen- 
tial bias the mean squared error grows exponentially with 
the dissipative work, it increases with Vw for the 1/P 
bias and linearly for thermodynamic integration. 

As can be seen in Fig. [2] in the slow switching regime 
the errors decrease with decreasing v in all four cases. To 
linear approximation in v, the error for straightforward 
fast switching, work biased umbrella sampling with the 
exponential bias 7r(W) = exp(— (3W/2), and work biased 
thermodynamic integration is identical and given by 



: N 



crjv = 2jk B TLv 
N N 



(81) 



For work biased umbrella sampling with bias n(W) = 
1/P(W) the error is smaller and the dependence on v 
stronger: 



/3a 3 



-JV 



{2 1 Lvf' 2 ^k B T 



(82) 



Here, the interval size U = W max — W m - m depends on the 
selection of the limits Wmax and W m i n . As a measure for 



As one proceeds to larger switching rates, the error grows 
most rapidly for straightforward fast switching. For work 
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FIG. 2: Error measure k 2 = e 2 N N for the particle pulled 
through a viscous fluid as a function of switching rate v ob- 
tained for straightforward fast switching (FS), work biased 
thermodynamic integration (TI), and work biased umbrella 
sampling with the umbrella functions ir = exp(— /3W/2) and 
7T = 1/P(W). Here, k B T = 1, 7 = 1, k = 1 and L = 5. 




V 



FIG. 3: Computational cost Ccpu = f?t for the particle 
pulled through a viscous fluid as a function of switching rate v 
obtained for straightforward fast switching (FS), work biased 
thermodynamic integration (TI), and work biased umbrella 
sampling with the umbrella functions ir = exp(— /3W/2) and 
ir = 1/P(W). Here, k B T = 1, 7 = 1, k = 1 and L = 5. 



biased umbrella sampling with the exponential umbrella 
function the growth is less pronounced and occurs at 
higher switching rates. For the two remaining cases (work 
biased umbrella sampling with tt = 1/P(W) and work 
biased thermodynamic integration) the growth of the er- 
ror slows down rather then speeding up with increasing 
switching rate. Since the work variance becomes con- 
stant for large v, the error reaches a constant value in 
the fast switching regime in all four cases. In this limit 
we obtain 



N 



k 2 T 2 9 9 
N 



for straightforward fast switching, 
9 L k 



N 



for work biased thermodynamic integration, 



N 



_ 2e P 2 L 2 k/i 



(83) 



(84) 



(85) 



for work biased umbrella sampling with the exponential 
umbrella function n(W) = exp(— (3W/2) and 



-N 



klT 2 yjJIFk 



N 



(86) 



in the case where the umbrella function ir(W) = 
1/P(W). Note, in Fig. [2j that these constant values 
may differ by orders of magnitude. 

The computational cost Ccpu determined for the four 
algorithms studied here is shown as a function of the 
switching rate in Fig. [3] In the slow switching regime the 
computational cost becomes constant for straightforward 
fast switching, work biased umbrella sampling with expo- 
nential bias and work biased thermodynamic integration. 



This behavior indicates that in these cases running few 
long trajectories or many shorter trajectories yields the 
same efficiency The situation is different for work biased 
umbrella sampling with the 1/P(W) umbrella function. 
In this case, Ccpu keeps decreasing for decreasing v even 
in the slow switching regime. For larger switching rates 
straightforward fast switching and work biased umbrella 
sampling with exponential bias become very inefficient. 
In the other two cases, increasing the switching rate en- 
hances the efficiency With the appropriate bias the opti- 
mum efficiency is obtained in the instantaneous switching 
limit, in which case the path based sampling methods re- 
duce to ordinary umbrella sampling and thermodynamic 
integration. Thus, conventional free energy calculation 
methods are superior to fast switching methods for the 
model studied in this section. Similar results have been 
obtained earlier for more complex models (ill . [l2| . 

In addition to the expected error also the bias 6jv is an 
important quantity determining the accuracy of a free 
energy calculation. According to Eq. (|3l?l) the bias can 
be calculated using the expressions in Eqs. (|35|) to (|40[) . 
For the pulled particle with Gaussian work distribution 
one obtains 



kjiT , rj2 2 



1). 



(87) 



in the case of straightforward fast switching. As can be 
seen by direct evaluation of Eq. (|32|). all other algorithms 
are bias free, i.e., bpf = 0, if the work distribution is 
Gaussian. 

It is interesting that in the case of work biased sam- 
pling with exponential umbrella function, the property 
of vanishing bias is only given for the particular form 
ir(W) = cxp(— PW/2). To be more specific, one can con- 
sider umbrella functions tt(W) = exp(— X(3W), where A 
is a real number between and 1. In this case, error and 
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bias are given by 



: -iY 



Ml! 

N 

,(l-\) 2 f] 2 v 2 w 



and 



b N 



k B T f 
2N 



,(l-A) 2 /3 2 



= A /3 (7y 



(88) 



(89) 



respectively. The error from the above equation is a 
minimum for A = 1/2. For this value of the parameter 
A the bias vanishes. In this sense, the umbrella function 
ir(W) = cxp(— (3W/2) is the optimum umbrella function 
among all exponential functions of the form exp(— X(3W). 
Note, however, that this property holds strictly only for 
processes with Gaussian work distribution. 

We now briefly return to the question of whether in 
expression (|48|) for the number of trajectories required to 
obtain an error of about k^T the average of the expo- 
nential can be safely replaced by the exponential of the 
average (see Eq. [50)) . For the process considered in this 
section, the dissipative work as well as the work distribu- 
tions are the same in the forward and reverse direction. 
Furthermore, W d = W because AF = 0. Hence, Eq. 
(1481) reduces to: 



N = exp(/3 2 a^) 



exp(2/3W), 



(90) 



where we have used Eq. (|73| which is valid only for Gaus- 
sian work distributions. Thus, in this case, the estimate 
N w exp(/3W) may be off by orders of magnitude. If, for 
instance, the correct estimate is N = 10.000, replacing 
the average of the exponential with the exponential of 
the average yields an estimate too low by a factor of 100. 



B. Expansion of ideal gas 

1. Model and work distribution 



limit of a piston that moves infinitely slowly. In this case, 
the distribution is a rapidly fluctuating function with an 
exponential envelope. In the fast switching limit, the 
work distribution acquires a growing delta-peak at zero 
work and small contributions at large work values. These 
large but rare work values are the dominant contributions 
to the exponential work average. 

In their work, Lua and Grosberg [lj| [34[ considered a 
single ideal gas particle in a cylinder with cross section 
area A sealed by a piston on one side. The volume avail- 
able to the gas particle is V = ALo, where Lq is the initial 
length of the cylinder. Starting from initial conditions 
distributed canonically with temperature T = (fcs/?) -1 , 
the volume of the cylinder is then increased by moving 
the piston from position Lq to a new position L in a 
time r with constant velocity v p = (L — Lq)/t. The 
ideal gas particle is assumed to collide elastically with 
the walls of the cylinder and the moving piston. At each 
clastic collision of the particle with the piston, a certain 
amount of kinetic energy is transferred from the parti- 
cle to the piston. This adiabatic case is in contrast to 
that considered in Ref . [28[ , where the system is coupled 
to a thermal bath and the velocities are taken from an 
equilibrium distribution at each collision with the piston. 
The sum of these energy transfers during one expansion 
is the work W carried out by the gas. Note that here, as 
in Lua and Grosberg's work, W denotes the work done 
by the system and not on the system. Accordingly, we let 
AF = F(A) — F(B) denote the free energy difference be- 
tween the initial and the final state. The maximum work 
theorem then implies that the average work performed 
by the system is bounded from above by the free energy 
difference, (W) < AF. With this sign convention, the 
Jarzynski equality is given by (exp(/3W)) — exp(fiAF). 

Lua and Grosberg have analyzed the statistics of the 
particle-piston collisions in detail for the special case of 
inverse temperature /3 = 1, mass of the particle m = 1 
and switching time r = 1. In this paper we keep the 
explicit dependence on all variables and obtain the fol- 
lowing work distribution : 



In this section, we will determine the efficiency of fast 
switching methods for the expansion of an ideal gas in a 
cylinder with moving piston. The work distribution for 
this process has been determined analytically by Lua and 
Grosberg [l|| [34| , who demonstrated that the Jarzynski 
identity holds exactly for arbitrary piston speed. In the 
limit of very rapid expansion this seems surprising as 
almost no gas particles collide with the piston performing 
work on it. Nevertheless, as shown by Lua and Grosberg, 
the tails of the Maxwcll-Boltzmann distribution provide 
a sufficient number of fast moving particles that perform 
exactly the amount of work required for the Jarzynski 
identity to be valid. 

The reason why we study the ideal gas expansion here 
is that its work distribution is non-Gaussian. As the sys- 
tem is not coupled to a heat bath during the transition, 
the work distribution remains non-Gaussian even in the 



P(W,L ,v p ,/3) = P S{W) + 



e 2 2m ™p' f(W,L ,Vp,P). (91) 

2ir mnVp 



where 




dve 



<(„- 



■ + l)v p 



(92) 

is the probability that the particle does not collide with 
the moving piston during the expansion process. This 
probability increases with increasing piston velocity v p 
leading to a growing delta-peak at W — in the work 
distribution. As shown by Lua and Grosberg, the number 
n of collisions between the moving ideal gas particle and 
the piston is completely determined by the work W done 
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by the gas during the entire expansion and is given by 
n(W,Lo,v p ) = 




2W 



(93) 



where [■ ■ 'J indicates the floor function that gives the 



largest integer less or equal to the argument. The func- 
tion f(W,Lo,v p ,P) appearing in Eq. (|9"Tj). called the 
overlap factor by Lua and Grosbcrg, is given by the piece- 
wise linear function 



f(W) = 



"(n -!)(! + ££) + 
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AT. (94) 



A variation of the piston velocity v p is equivalent to a 
variation of the temperature in the sense of the following 
scaling relation: 

P(W,L ,^,f3) = X 2 P(X 2 W,L Q ,v p ^). (95) 

This equation, which can be easily derived directly from 
Eq. (j9l"j) , holds, because due to the elastic hard collisions 
the same sequence of events occurs if both v and v p are 
scaled by the same factor. 
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FIG. 4: Average work (W) (solid red line) performed by the 
ideal gas during the expansion as a function of the piston 
velocity v v for Lo = 1, AL = 1 and j3 — 1. Also shown is 
the work done on the ideal gas during the compression (solid 
green line). The free energy difference AF — fceTlog[(Lo + 
AL)/Lo] — log(2) is shown as a dashed line. For v p — > 0, the 
average work converges towards (1/2)[1 - L§/(L + AL) 2 ] = 
0.375. 



The average work 

<W0 



P(W)WdW 



(96) 



performed during the expansion of the cylinder from 
Lo to L = Lq + AL calculated from the work distri- 
bution Eq. (|91[) by numerical integration is shown in 



Fig. 2J It is clear from the picture that the aver- 
age work (W) is lower than the free energy difference 
AL = F(L )~F(L) = /c B Tlog[(Lo+AL)/L ] even in the 
limit of infinitely slow switching. In this limit, the pro- 
cess studied by Lua and Grosberg can be viewed as the 
quasistatic adiabatic expansion of an ideal gas where the 
average work follows from basic thermodynamic equa- 
tions for the ideal gas: 



(wo 



PdL 



1 



-dL 



k B T 



1 



7 



1-7 



(97) 



where 7 = C v jCy = (/ + 2)// and / is the number of 
degrees of freedom. In our one-dimensional case, / = 1 
and 7 = 3 such that 



(W) 



k B T (L 



T 2 



L 2 



(98) 



For fast switching, on the other hand, the average work 
converges to zero because less and less collisions of the 
moving particle with the piston occur as the piston ve- 
locity is increased. Nevertheless, the Jarzynski equality 
holds over the whole range of piston velocities as can be 
directly verified from Eq. (|9ip . 

In the slow switching limit, i.e., for v p <C 
\J (k-BT/m)AL/(Lo + AL), the work distribution can 
also be obtained by considering adiabatic invariants 
[29l l30l . I3~i1 . [32l . [33( 1 . For a one-dimensional mechanical 
system the principle of adiabatic invariance states that 
the action variable 



(99) 



is conserved under slow variation of an external param- 
eter in the Hamiltonian. In the above equation the con- 
tour integral extends over one period of the motion for 
constant energy E = p 2 /2m. For the Lua-Grosberg 
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model this relation implies that the quantity EL 2 is con- 
served during a slow, i.e., adiabatic, transformation of 
the Hamiltonian. As a consequence, the energy E of the 
system at the end of an expansion from Lq to L is a 
unique function of the initial energy Eg: 



Eq. 



(100) 



Thus, the work carried out by the gas starting from any 
initial condition with energy Eq is given by W — Eq — 
(Lo/L) 2 E = E [1- (Lq/L) 2 }. The work distribution in 
the slow switching limit can then be obtained as average 
over the canonically distributed initial conditions: 



P(W) = 



ft 



2-Km 



dp e -P p2 / 2m 5 



2m \ L 2 



Integration yields 



P(W) 



/3L 2 
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tt{L 2 
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(101) 



(102) 



which is the one-particle one-dimensional special case of 
the work distribution derived by Crooks and Jarzynski 
for an ideal gas of TV particles [33]. It is easy to verify 
directly that for this work distribution that the Jarzynski 
equality holds and the average work (W) performed by 
the gas is given by Eq. ([98)) . Note that even tough the 
expansion occurs infinitely slowly in this case, the work 
distribution is not Gaussian and the width aw and the 
average (W) of the work distribution are related by 



aw 



k B T 
V2 



L 2 -L 2 



L 2 



= V2(W) 



(103) 



rather than by Eq. (|73|) . 
In the fast switching 
AL), 



limit, i.e. 



for 



» 



a/ (fceT '/ 'to) AL I (Lq + AL), only trajectories with or 
1 collision yield important contributions to the work dis- 
tribution. In this case, one obtains 13411: 



P(W) 



where Pq is given by 




(2mv';+W) 



Av^m 2 



Pn 



1 



1 



V27TTO/3 



AL 
1~o~ 



(104) 



(105) 



In the following paragraphs we will also consider the 
statistical errors for the reverse process, namely the com- 
pression of the ideal gas from the cylinder length Lq + AL 
to Lq. The work distribution Pr(W) for the reverse pro- 
cess is simply given by the Crooks relation Eq. (|4"Tj) . 
which, with the sign conventions for W and AF used in 
this section reads: 



Pr(-W) = P(W)e 



P(W-AF) 



P(W)e 



pw 



Lq 



(106) 



2. Efficiency 

We next determine the efficiency of fast switching 
methods for the ideal gas expansion from Lq to Lq + AL 
as well as for the compression from Lq + AL to Lq. In 
this case, we investigate the straight forward approach, 
biased sampling with cxp(/3W/2) umbrella function, and 
thermodynamic integration. We do not, however, con- 
sider the \jP bias here, because the work distribution 
for the ideal gas expansion includes a (5-peak that creates 
problems in this case. 

When estimating statistical errors of fast switching 
simulations another difficulty arises in some situations. 
As detailed in Sec. IIII| the mean square error e 2 N as well 
as the bias bjq can be calculated from integrals over the 
work distribution. In the case of the ideal gas expan- 
sion, however, some of the integrals from Eqs. (|35p - (|4"0|) 
diverge for certain expansion ratios. To see this, con- 
sider the mean square error estimated for a straightfor- 
ward fast switching simulation using the sign convention 
of this section: 



kjT 2 
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f e -2/3AF l^pWx 
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(107) 



The average (exp(2/?W)) = / dWP{W) exp(2/?IF) exists 
only if the work distribution P(W) decays sufficiently 
rapidly with increasing W to compensate for the fast 
growth of exp(2(3W). For large work values W, the work 
distribution for the ideal gas expansion behaves as 



P(W) ~ exp 



-pw 



Lr 



L 2 -L 2 



(108) 



where we have neglected all factors that are immaterial 
for the convergence of the integral. Thus, (exp(2f3W)) is 
finite only if 



L 2 



L 2 



T 2 



-2 > 0, 



(109) 



or, equivalently, if the ideal gas is expanded by no more 
than a factor of \2, (L/Lq) < \/2. If this condition, 
which is independent of the switching rate, is not met, 
X = cxjp([3W) has infinite variance and the expressions 
from Sec. IIIII can not be applied to estimate the statis- 
tical error of the free energy estimate. In this 
generalization of the law of large numbers [35| still guar- 
antees that the free energy converges to the correct value, 
but the fat-tailed work distribution prevents Gaussian 
statistics from arising even for large sample sizes. 

For umbrella sampling with the exponential bias 
cxp(/3Vy/2) a similar analysis of the average (Y 2 )^ re- 
veals that the corresponding integral converges only pro- 
vided that [L/ Lq) < vo. For larger expansion ratios, the 
integral diverges for any switching rate. In contrast, the 
integrals associated with the thermodynamic integration 
scheme always converge regardless of the expansion ratio 
L/Lq and the switching rate. 
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The situation is slightly different for the ideal gas com- 
pression. In this case, straightforward fast switching 
yields error estimates that are finite for any expansion 
ratio. Since for umbrella sampling with exponential bias 
cxp(pW/2) and thermodynamic integration the errors in 
forward and backward direction are identical, the same 
limitations as for the ideal gas expansion apply in these 
cases. 

In the following we will study the efficiency for ex- 
pansion ratios of L/Lq = 1.2 and L/Lq = 2.0. For 
L/Lq = 1.2 the error estimates for all fast switching 
methods are finite both for the expansion and the com- 
pression. For L/Lq = 2 we consider only the ther- 
modynamic integration method and straightforward fast 
switching for the ideal gas compression. All other meth- 
ods yield diverging error estimates for this expansion ra- 
tio. 



1 
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— exp(W/2) 
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FIG. 5: Error measure k 2 = e%N as a function of the piston 
velocity v p for different fast switching methods and an expan- 
sion ratio of L/Lo = 1.2. The black lines denote the straight 
forward fast switching errors for expansion (solid line) and 
compression (dashed line). 

Figure [5] shows the errors as a function of the piston ve- 
locity v p obtained by numerical integration of Eqs. (|35|) - 
PP)) for an expansion from Lq = 1 to L = Lq + AL = 1.2. 
Here and in the subsequent calculations the inverse tem- 
perature of the initial state was (3 = 1. For infinitely slow 
and infinitely fast switching the errors can be calculated 
analytically and were found to agree perfectly with the 
results obtain by numerical integration. As can be seen 
in Fig. the errors become independent of the switching 
rate v p for all three methods in the slow switching limit. 
In the fast switching limit, the error grows rapidly with v p 
particularly for straightforward fast switching and the ex- 
ponential bias. The thermodynamic integration scheme, 
for which the growth is slower, yields smaller errors for 
all switching rates. For all three cases the errors diverge 
for v p — > oo. The overall smallest errors, however, are ob- 
tained if the reverse process, i.e., the compression, is sim- 
ulated with straightforward fast switching (dashed line in 
Fig. O. In the fast switching limit this error becomes 



AL/Lq. Since the initial position of the particle is uni- 
formly distributed in the interval [0, L], this value is sim- 
ply the probability that the fast moving piston hits the 
particle at least once. Note that for thermodynamic inte- 
gration and umbrella sampling with the exponential bias 
the errors arc identical for expansion and compression. 

From the errors depicted in Fig. one can determine 
the computational cost required to obtain an accuracy 
of ksT in the free energy estimate [see Eq. (pT)) ]. This 
computational cost is shown in Fig. [6] for different fast 
switching methods and an expansion ratio of L/Lq = 1.2 
as a function of the switching rate. For straightforward 
fast switching in the expansion direction, thermodynamic 
integration and umbrella sampling with the exponential 
bias an optimum switching rate of about v p ~ 1 exists. 
In the case of straightforward fast switching in the com- 
pression direction, however, no such minimum exists and 
the computational cost keeps decreasing for increasing 
switching rate. 




io" 2 10" 1 10° io 1 

v 

p 

FIG. 6: Computational cost Ccpu = k 2 t as a function of 
the switching rate v p for different fast switching methods and 
and expansion ratio of L/Lq = 1.2. The black lines denote 
the straight forward fast switching errors for expansion (solid 
line) and compression (dashed line). 

As explained above, for an expansion ratio of L/Lo = 2 
our expressions yield finite statistical errors only for the 
thermodynamic integration method and for straightfor- 
ward fast switching in the compression direction. The er- 
rors for these cases are shown in Fig. [7J Straightforward 
fast switching of the ideal gas compression yields errors 
that are smaller than the thermodynamic integration er- 
rors for all piston velocities. Both in the slow and fast 
switching limit the straightforward fast switching errors 
become constant. Accordingly, the computational cost 
of the straightforward compression is decreasing as 1 /v p 
in the fast switching limit while the thermodynamic in- 
tegration error goes through a minimum at about v p ~ 1 
and then rapidly grows for increasing piston velocity v p 
(see Fig. EJ). Thus both for L/Lq = 1.2 and L/Lq = 2.0 
the straightforward fast switching simulation of the ideal 
gas compression is the most efficient among the methods 
studied here. However, optimum efficiency is obtained 
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for infinitely fast compression, in which case straightfor- 
ward fast switching reduces to Zwanzig's perturbative 
approach applied to configurations rather than trajecto- 
ries GB. 




FIG. 7: Error measure e N N as a function of the piston ve- 
locity v p for the ideal gas expansion/compression with an ex- 
pansion ration of L/Lo = 2.0. The dashed black line denotes 
the straightforward fast switching errors (FS) and the solid 
red line the thermodynamic integration errors (TI). 




FIG. 8: Computational cost Ccpu as a function of the switch- 
ing rate v p for the ideal gas expansion/compression with an 
expansion ratio of L/Lo = 2. The dashed black line denotes 
the straight forward fast switching results (FS) for the com- 
pression and the solid red line the thermodynamic integration 
results (TI). 



3. Many-particle gas 

The model discussed in the previous section consisted 
of a single particle in one dimension. In this section, 
we study the expansion of an ideal gas of M particles 
and analyze how the error depends on the number of 
particles. The system of strictly non-interacting particles 
considered here is similar, but not identical to the dilute 



gas of weakly interacting particles studied by Jarzynski 
and Crooks in Ref. [36] ]. For the M-particle ideal gas, 
the work performed on the piston during the expansion 
is the sum of the contributions of each particle. 



W = ^Wi. 



(110) 



i=l 



Here Wi is the work performed on the system by particle 
i. According to Eq. (|34|) . the statistical error in the 
free energy calculated with straightforward fast switching 
from N repetitions of the process is given by 



-N 



N 



(xy 



- 1 



(in) 



For the M-particle system, the averages (X) and (X" 2 
are 



(X) 



and 



(X 2 ) = ( e -^i 2Wi ) = (X 2 ) 



2\M 



(112) 



(113) 



where (X\) and (X 2 ) are the corresponding averages of 
the one-particle system. Here, we have used the fact that 
work contributions of the single particles are identical 
and statistically independent. The error €n(M) for the 
A/-particle system becomes 



■ N 



(M) 



k 2 B T 2 
N 



(Xi) 2 



M 



(114) 



The fraction (Xf) j \X\) 2 is always larger than 1 and 
therefore the error e diverges exponentially with the par- 
ticle number in the case of straightforward fast switching. 

A similar analysis can be carried out for the error of 
work biased umbrella sampling simulations. For an um- 
brella function 7r(VF) that factorizes into a product of 
single particle umbrella functions 7Ti(Wj) (this is for in- 
stance the case for the exponential bias exp(— (3W/2)), 



the error is given by: 
4(M) 



(115) 



k 2 B T 2 
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(Xi)l 
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V (Yi)l 



M 



(X^ 



M 



(116) 



Here, as in the previous paragraph, the averages are 
single-particle averages. For large M, this expres- 
sion will be dominated by either ((X 2 )^/ (Xi) 2 T ) M or 
(P?)*/(ri)D M depending on whether {Xf) v /{X 1 )l or 
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(Yi) 7T /(Yi)^ r is larger. (It is easy to see that the larger 
one of (Xi) 7r /(X 1 )l and (Yj 2 ) n / '(Yi) 2 - is also larger than 
(XiYi) 7r /(Xi) w (Y'i) 7r , which is always positive.) Thus, 
also in this case, the error in the free energy grows expo- 
nentially with the system size M. 

For the work biased thermodynamic integration the er- 
ror results from an integral over the work variance at dif- 
ferent values of an auxiliary parameter a (see Eq. (|57t ). 
Since in the M-particle ideal gas the work W is a sum of 
M independent contributions, the work variance of the 
entire system is M times the single particle work vari- 
ance for each value of a. Hence, the M-particle error 
e(M) scales with the square root of the number of parti- 
cles, 

e%(M) = Me 2 N (l), (117) 

where £at(1) is the error for the one-particle gas. Due 
to this kind of scaling, the thermodynamic integration 
method is superior to all other fast switching methods 
studied here for a strictly non-interacting ideal gas. 

VI. CONCLUSION 

For irreversible transformations, the exponential aver- 
age of Eq. ([2]) is dominated by rare, but important work 
values causing possibly large statistical inaccuracies in 
fast switching free energy computations. Several path 
sampling methods have proposed that aim at reducing 
these inaccuracies by enhanced sampling of those trajec- 
tories that contribute most to the average. In this paper, 
we investigate how these methods perform in the case 
of two simple models for which the work distribution is 
known analytically. 

For Gaussian work distributions, the computational 
effort required to achieve a given accuracy in the free 
energy can be determined analytically. In the case of 
a particle pulled through a viscous fluid, flat histogram 
sampling as well as Sun's work biased thermodynamic 
integration perform best. For both methods optimum 
efficiency is obtained in the limit of infinitely fast switch- 
ing in which case these path sampling methods reduce 
to conventional flat histogram sampling and thermody- 
namic integration in configuration space rather than in 
trajectory space. 

The errors obtained in the case of the ideal gas expan- 
sion/compression, for which the work distribution was 
determined analytically by Lua and Grosberg [l6|, [34| . 
show a qualitatively different behavior. In contrast to 
the particle pulled through a fluid, the work distribution 
for the expansion/compression process does not approach 
a constant shape in the fast switching limit. Rather, the 
work distribution keeps changing even for very large pis- 
ton velocities. As a consequence, the computational cost 
of all methods goes through a minimum for intermedi- 
ate piston velocities and grows for rapid expansion. The 
only exception here is the straightforward fast switching 
simulation of the gas compression, which is superior to 



the path sampling methods for all piston speeds. In this 
case, the computational cost keeps decreasing even for 
very large piston speeds where the fast switching simu- 
lation reduces to Zwanzig's pcrturbative approach. It is 
interesting that for certain expansion ratios only thermo- 
dynamic integration and straightforward fast switching of 
the compression process yield finite error estimates, while 
all other method yield diverging errors. This behavior is 
due to the fat tails of the work distribution that occur 
for larger expansion ratios. 

The superiority of straightforward fast switching for 
the ideal gas compression is most likely peculiar to the 
one-particle case. In the many-particle gas, the errors 
scale exponentially with particle number for all methods 
except for thermodynamic integration, in which case they 
scale with the square root of the particle number. Thus, 
for an ideal gas consisting of many particles Sun's work 
biased thermodynamic integration method seems to of- 
fer the most efficient free energy calculations. Since for 
this method the efficiency is maximum at intermediate 
piston velocities, this is an example where a fast switch- 
ing method for the calculation of free energies is more 
efficient than conventional configuration space methods. 

In general, the efficiency of a fast switching simulation 
depends on whether the process is carried out in forward 
or reverse direction. In particular, a straightforward fast 
switching simulation is more efficient in the direction in 
which more work is dissipated [l8j as confirmed here for 
the Lua-Grosberg ideal gas expansion and compression. 
For all the path sampling methods studied here, how- 
ever, it can be shown using the Crooks identity that the 
work distribution is effectively symmetrized in a way that 
makes the corresponding efficiency identical in both di- 
rections. 

In summary, path sampling methods such as work bi- 
ased thermodynamic integration or work biased umbrella 
sampling may help to reduce the statistical errors en- 
countered in the calculation of equilibrium free energies 
from non-equilibrium simulations. The examples stud- 
ied in this paper indicate that often conventional free 
energy estimation methods are likely to be superior to 
fast switching simulation. Nevertheless, advanced fast 
switching methods such as Sun's work biased thermo- 
dynamic integration may be competitive particularly if 
a large number of degrees of freedom significantly con- 
tribute to the free energy change that one wants to com- 
pute. 
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